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Introduction. 

The  ReaxFF  2008  C/H/0  force  field1,  popularly  known  as  the  Chenoweth  force  field 
has  been  successfully  applied  to  study  the  oxidation  of  high-molecular  weight  fuels. 
However,  while  this  force  field  has  been  successfully  used  for  studying  the  initial 
stages  of  pyrolysis  and  combustion  for  complex  fuels  [e.g. 210)  the  performance  of 
this  force  field  in  studying  smaller  hydrocarbon  combustion  is  not  satisfactory, 
especially  for  reactions  involving  carbon  monoxide.  The  main  goal  of  our  project  is 
to  improve  the  Chenoweth  force  field  for  Cl  chemistry  so  that  we  can  capture  all  the 
important  species  mentioned  in  literature  [Li  et  al,  J.  of  Chem.  Kinet,  2008 )  and 
along  with  the  reactions  involving  those  species,  while  retaining  the  quality  of  the 
force  field  for  larger  hydrocarbon  pyrolysis  and  combustion  reactions. 

Furthermore,  in  order  to  extend  communication  between  experiment  and 
simulations,  we  have  developed  a  computational  strategy  for  evaluating  flame  speed 
in  fuel/oxidant  mixtures. 


Method  development. 

In  order  to  improve  ReaxFF  description  for  Cl-chemistry  -  more  specifically,  the 
CO-chemistry,  we  modified  the  ReaxFF  angular  term  that  relates  to  over¬ 
coordination  effects  on  the  central  atom. 
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Equation  1:  Modified  ReaxFF  angular  function  for  central-atom  overcoordination. 


H-C-O  Angle  (Degree) 


Figure  1:  Comparison  of  the  Chenoweth  and  the  C/H/O  2015  ReaxFF  description  of 
H-C=0  angle  bending  in  a  HCO-radical  against  quantum  mechanical  [QM)-data.  QM 
calculations  were  performed  at  the  DFT  6-311G**/B3LYP  level. 
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As  Figure  1  indicates  -  this  modification  enables  our  new  ReaxFF  C/H/O  description 
(C/H/O  2015]  to  properly  capture  the  H-C=0  equilibrium  angle  in  the  HCO-radical  - 
which  is  a  key  intermediate  in  Cl  chemistry.  In  the  Chenoweth  description,  this 
angle  is  forced  to  be  linear  [Figure  1],  but,  according  to  quantum  mechanical  [QM) 
calculations,  it  should  have  an  equilibrium  angle  of  around  120  degrees.  As  Figure  1 
indicates,  the  modification  to  the  ReaxFF  angular  function  enables  C/H/O  2015  to 
properly  reproduce  the  C/H/O  equilibrium  angle  and  is  distortion  energies. 


Training  of  the  C/H/O  2015  parameters. 

After  introducing  the  new  angular  term,  we  re-trained  the  ReaxFF  C/O/H  bond, 
angular  and  torsion  terms  against  the  HCO-chemistry,  as  revealed  from  QM  [Density 
Functional  Theory  [DFT)-calculations  by  Hsu  and  co-workers11.  Figure  2  compares 
the  DFT  data  from  Hsu  et  al.  to  the  ReaxFF  Chenoweth  and  C/H/O  2015  [Dec.14) 
results  -  showing  that  the  new  ReaxFF  parameters  provide  a  substantially  improved 
description  of  the  HCO-related  chemistry. 


Figure  2:  Comparison  of  the  Chenoweth  et  al.  and  the  C/H/O  2015  [Dec.  2014) 
ReaxFF  descriptions  of  the  HCO  +  O2  reaction  to  CO  +  HO2  with  the  DFT-data  from 
Hsu  et  al.  [1996). 


We  also  included  other  HCO-reaction  channels,  including  HCO+O2  OH  +  CO2  and 
alternative  HCO+O2  CO  +  HO2  pathways  in  our  training  and  found  that  the  C/H/O 
2015  parameters  are  in  overall  good  agreement  with  DFT  for  all  these  reactions. 
Subsequently,  we  added  the  other  ReaxFF  training  items  using  in  the  Chenoweth  et 
al.  development  to  the  HCO-training  set  -  thus  ensuring  that  the  C/H/O  2015  force 
field  retains  the  Chenoweth  et  al.  quality  for  more  complex  fuels. 


Validation  of  the  C/H/O  2015  parameters  using  high-temperature  molecular 
dynamics  simulations. 

In  order  to  validate  the  new  C/H/O  2015  ReaxFF  description  we  performed  a  series 
of  high-temperature  [3000K)  reactive  molecular  dynamics  simulations  on  H2/O2/CO 
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mixtures.  From  these  simulations  we  could  directly  evaluate  the  available  reaction 
channels  for  conversion  of  CO  to  CO2.  Figure  3  shows  the  reaction  channels,  as 
observed  from  these  simulations,  using  the  Chenoweth  and  the  C/H/O  2015  force 
fields: 


Figure  3:  CO  to  CO2  reaction  channels  in  a  CO/H2/O2  mixture,  as  observed  from 
ReaxFF  simulations  using  the  Chenoweth  [top)  and  the  C/H/O  2015  (bottom)  force 
fields. 

The  Chenoweth  force  field  only  captures  a  single  channel  for  CO-oxidation  -  by 
direct  oxidation  using  either  atomic  oxygen,  HO,  or  HO2  as  an  oxidant.  However,  the 
C/H/O  2015  ReaxFF  description  recognizes  additional  channels  through  an  HCO 
intermediate.  Table  1  summarizes  the  species  observed  in  the  Chenoweth  and 
C/H/0/2015  simulations,  and  compares  them  against  the  available  experimental 
data,  clearly  demonstrating  that  the  C/H/O  2015  gives  a  far  more  complete 
description  of  the  CO-reaction  channels. 
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Table  1:  Comparison  of  reactions  observed  during  the  ReaxFF  simulations  using  the 
Chenoweth  and  the  C/H/O  2015  force  field  and  experimental  data  (1:  Li  et  al., 
J.Chem.Kinet.  39:3,  2007,  2:  Mueller  et  al.  J.Chem.Kinet.  31:113, 1999,  3:  Starik  et  al., 
Combustion,  Explosion  and  ShockWaves,  46:5,  2010) 

Atomistic-scale  simulations  of  flame  speed. 

In  order  to  develop  new  concepts  for  comparing  atomistic-scale  simulations  with 
experiment  we  have  designed  a  simulation  strategy  for  measuring  flame  speeds 
from  ReaxFF-based  molecular  dynamics  simulations.  For  now,  we  have  been  using 
the  Chenoweth  et  al.  ReaxFF  description  in  these  simulations  -  we  are  planning  to 
transition  to  the  C/H/O  2015  description  in  the  near  future. 

Initially,  we  set  up  the  simulations  in  a  relatively  small  system  (Figure  4),  where  we 
equilibrate  a  02/fuel  mixture  (here,  the  fuel  is  n-butane)  to  a  temperature  of  1800K. 
At  this  temperature,  on  the  ReaxFF  pico-to  nanosecond  time  scale,  this  mixture  is 
non-reactive.  Subsequently,  we  replace  2  O2  molecules  with  O-radicals  and  perform 
a  NVE  (microcanonical)  ensemble  simulations.  The  O-radicals  initiate  a  combustion 
reaction  that  spreads  through  the  simulation  box  -  raising  temperature  (Figure  5) 
and  causing  radical  species  to  spread  through  the  system.  By  monitoring  the  OH- 
radical  front  as  a  function  of  time,  we  can  derive  a  flame  speed  from  our  simulations. 
As  indicated  in  Figure  5  -  we  successfully  identified  first-stage  combustion  - 
primarily  resulting  in  formaldehyde  formation  -  and  second-stage  combustion  - 
resulting  in  complete  oxidation  to  CO2  and  H2O  -  during  our  simulations,  where  the 
first-stage  is  associated  with  a  limited  temperature  increase,  following  by  a  much 
more  substantial  temperature  increase  during  the  second  stage. 
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Figure  4.  Mixture  of  12  Butane/78  O2  in  a  12x12x80  Angstrom  periodic  box,  as  used 
in  the  small-scale  flame-speed  ReaxFF  simulations.  The  O-atoms,  intiating  the 
combustion,  were  added  on  the  extreme  left  side  of  the  box. 


Time  (ns) 


Figure  5:  Main  species  and  system  temperature  observed  during  the  small-scale 
(Figure  4)  ReaxFF  NVE-simulations  on  combustion  initiations.  We  observe  a  distinct 
maximum  in  formaldehyde  concentration,  connected  to  a  modest  temperature 
increase,  around  t=0.6  ns.,  followed  by  a  complete  combustion  and  a  far  more 
significant  temperature  rise  during  the  later  stages  of  the  simulation. 


Our  analyses  from  the  small-scale  simulations  (Figure  4)  showed  a  very  fast  OH- 
radical  diffusion,  making  flame  speed  measurements  highly  uncertain.  For  that 
reason,  we  increased  the  box  by  a  factor  25  in  the  z-direction  (12x12x2000 
Angstrom)  allowing  more  straightforward  observation  of  the  OH-radical  profiles  as 
a  function  of  time. 
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Figure  6  shows  the  OH-concentration  profiles  as  a  function  of  time  for  the  2000 
Angstrom  box  -  showing  a  distinct  OH-profile,  from  which  a  flame  speed  can  be 
derived. 
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Figure  6:  OH-concentration  profiles  as  a  function  of  time  for  the  12x12x2000 
Angstrom  box  -  all  orange  dots  indicate  the  presence  of  an  OH-radical.  Note  that, 
due  to  the  periodicity,  the  OH-front  advances  on  both  sides  of  the  periodic  box. 


From  the  OH-profiles  we  could  estimate  a  flame  speed  of  around  125  m/s.  This  is 
significantly  faster  than  the  experimental  flame  speeds  for  similar 
hydrocarbon/oxidizer  mixtures  (0.4  m/s),  however,  both  the  ReaxF  initial 
temperatures  (1800K)  and  pressures  (55MPa)  are  substantially  higher  that 
experiment  (T=300K,  P=latm).  Larger-scale  and  longer-timescale  ReaxFF 
simulations  will  allow  us  to  lower  the  initial  pressure  and  get  closer  to  the 
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experimental  conditions.  However,  even  with  the  differences  between  simulation 
and  experiment,  we  can  already  use  this  tool  to  evaluate  fuel/oxidizer  ratio  and  fuel 
structure  impact  on  flame  speed  -  and  connect  this  with  experimental  trends. 
Furthermore,  the  high-pressure  ReaxFF  simulations  may  provide  important 
information  for  future  experiments  on  high-pressure  combustion. 

H-abstraction  by  O2  from  Alkane 

H-abstraction  from  alkanes  by  O2  molecule  is  a  very  important  reaction  in 
combustion,  especially  for  smaller  hydrocarbons.  Smaller  hydrocarbons  do  not 
pyrolyze,  therefore  H-abstraction  by  a  nearby  O2  molecule  is  the  most  probable 
initiation  reaction.  It  has  been  previously  reported  that  the  reaction  CH4  +  02=  CH3 
+  HO2  has  a  reaction  barrier  of  ~50  kcal/mol  (Fig  1(a))  in  the  C/H/0  2008 
description.  The  experimental  barrier  for  this  reaction  is  52  kcal/mol12.  However,  it 
has  been  observed  that  this  reaction  takes  place  very  easily  while  simulating 
methane/02  system  at  2500K  in  spite  of  having  a  very  high  reaction  barrier.  From 
the  simulation  trajectory,  a  transition  state  for  this  reaction  is  captured  (Fig  1(b)), 
which  gives  a  reaction  barrier  of  only  ~31  kcal/mol.  This  structure  was  not  possible 
to  capture  through  the  bond  restraint  method  as  both  the  0-atoms  of  the  O2 
molecule  makes  weak  bonds  with  the  nearby  hydrogen,  which  is  about  to  be 
abstracted. 

However,  our  newly  developed  force  field  C/H/0  2016  estimates  an  energy  barrier 
of  ~48  kcal/mol  (TS  structure  shown  in  Fig  1(c)).  In  order  to  verify  the 
nonexistence  of  any  other  lower  energy  transition  state  structure  -  we  ran  sample 
simulations  of  the  same  methane/02  system  at  2500K.  During  the  simulation,  lesser 
amount  of  H-abstraction  by  O2  is  observed.  Thus,  now  we  have  to  simulate  our 
system  at  a  higher  temperature  (2500-3000K)  to  capture  methane  chemistry  at  a 
shorter  time  scale,  as  the  initiation  is  much  slower  now.  Previous  ReaxFF  study 
shows  that,  after  this  initiation  reaction,  the  two  products  of  this  reaction  recombine 
and  generates  CH3OOH,  which  then  produces  CH2O.  The  recombination  of  CH3  and 
HO2  is  only  possible  at  lower  temperature  as  these  two  will  stay  close  to  one 
another.  However,  C/H/O  2016  gives  a  slightly  different  pathway  of  CH2O 
generation  at  higher  temperature.  At  that  temperature,  OH  and  0  radicals  are 
generated  by  the  0-0  bond  dissociation  of  HO2.  This  0  radical  attacks  the  methyl 
radical  and  produces  CH3O,  which  then  leads  to  CH2O  formation.  Thus,  a  different 
pathway  of  CH2O  from  methane  combustion  is  captured  using  C/H/0  2016  force 
field. 
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Number  of  Molecules 


Fig  1:  The  transition  state  structure  of  the  reaction  CFU  +  02=  CH3  +  HO2  using 
C/H/0  2008  description  (a)  through  bond  restrain  (b]  through  simulation,  and  (c) 
through  bond  restrain  using  C/H/0  2016  force  field 

Product  Distribution  from  n-Butane  oxidation 
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Figure  2:  Time  distribution  of  temperature  and  mass  fraction  of  fuel  (C4H10], 
important  intermediate  species  (OH  and  CH2O]  and  combustion  products  (CO2  and 
H20]  using  (a)  C/H/0  2008  and  (b)  C/H/0  2016  force  field. 


Previously,  we  mentioned  that  we  have  developed  a  mechanism  to  quantify  laminar 
flame  speed  related  quantity  like  ignition  front  speed  using  ReaxFF.  At  that  time,  we 
used  C/H/0  2008  description  in  the  simulations  with  a  future  goal  to  transfer  to  the 
newly  developed  version.  Figure  2  shows  a  comparison  of  temperature  and  species 
distribution  from  simulations  using  C/H/0  2008  and  2016  description  -  same 
procedure  was  followed  for  both  the  cases.  We  initially  put  12  n-Butane  and  78  O2 
molecules  in  a  box  of  12A  x  12  A  x  80  A.  The  system  is  them  minimized,  equilibrated 
at  1800K  using  NVT-MD.  Next  NVE-MD  simulation  is  performed  at  1800K, which  will 
allow  us  to  observe  a  temperature  jump  indicating  combustion.  The  system  is 
initialized  with  two  0  atoms  to  bypass  the  high  ignition  delay  time.  Figure  2  shows 
that  the  temperature  and  species  distribution  profile  for  both  the  description  are 
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similar,  except  the  distribution  of  formaldehyde  (CH2O).  In  C/H/O  2008  force  field,  a 
clear  formaldehyde  peak  is  visible  before  the  OH  radical  build  up.  We  know  that  this 
is  a  characteristic  of  low  temperature  combustion  and  should  not  happen  at  our 
simulation  temperature  [1800K).  However,  in  C/H/O  2016  description,  no  clear 
peak  of  CH2O  is  observed  before  the  ignition  takes  places,  also  the  maximum  amount 
of  CH2O  is  almost  50%  of  the  C/H/O  2008  case.  As  discussed  previously,  C/H/O 
2008  has  much  lower  reaction  barrier  of  H-abstraction  by  O2,  which  is  one  of  the 
primary  paths  of  CH2O  formation.  This  problem  is  now  addressed  in  the  new  force 
field,  which  gives  a  lower  CH2O  yield.  Also,  no  relation  is  observed  between  CH2O 
production  and  system  temperature  indicating  no  effect  of  low  temperature 
chemistry. 


Initiation  Mechanism  for  Pyrolysis  and  Oxidation  of  Large  Hydrocarbons 

In  order  to  verify  our  newly  develop  force  field  is  performing  well  for  large 
hydrocarbons;  we  run  some  simulations  to  study  the  pyrolysis  and  oxidation  of  n- 
pentane.  To  get  some  statistically  meaningful  results,  we  perform  10  independent 
simulations  of  n-pentane/Ch  system.  At  first,  we  put  10  n-pentane  and  80  O2 
molecules  in  a  cubic  box  of  25A  side  length.  The  simulation  box  is  first  energy 
minimized  and  then  equilibrated  at  1500K  with  reactions  between  C-0  and  H-0 
turned  off  in  order  to  prevent  any  reaction.  Then  10  different  starting  points  are 
chosen  and  simulations  are  run  at  2500K  with  full  reactivity. 

In  all  the  10  simulations  we  perform,  the  system  initiated  through  pyrolysis  only, 
not  a  single  initiation  reaction  involving  H-abstraction  by  O2  is  observed.  This 
observation  is  consistent  with  our  previous  discussion  as  now  the  barrier  for  H- 
abstraction  reaction  by  O2  has  a  much  higher  value,  close  to  the  experimental  one. 
The  initial  decomposition  of  n-pentane  took  place  mainly  in  the  following  two 
pathways: 

H3°^^^^^^/CH3  - ^  ^^ch2  cr^3 


■>> 


CH„ 


Table  1  lists  the  percentage  of  each  reaction  observed.  Reaction  (1]  dominates  over 
reaction  [2] -this  trend  also  supports  a  previous  study10. 


Number 

Products 

%  of  Reaction 

ReaxFF  C/H/O  2016 

Allara  et  al.13 

1 

\^H2 

■^H3 

H2 

60 

62.6 

2 

CH3 

40 

22.3 

Table  1:  Percentage  of  the  different  reactions  observed 


Similarly,  we  performed  10  independent  simulations  with  alkyl  benzene  [butyl 
benzene)  to  see  how  the  new  force  field  performs.  Experimentally  it  is  observed  that 
pyrolysis  of  alkyl  benzenes  initiate  trough  beta  cleavage,  so  we  expect  to  see  high 
percentage  of  the  following  reaction 


ch3 


rrCH‘ 
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In  stead  we  are  observing  a  high  percentage  [50%)  of  the  alpha  cleavage  breaking  in 
our  simulation.  - ^ 


Thus,  the  newly  developed  force  field  has  failed  to  capture  the  kinetics  of  alkyl 
benzene  pyrolysis;  further  training  is  required  to  address  this  issue. 


Conclusions. 

We  have  successfully  improved  the  2008  ReaxFF  Chenoweth-combustion 
description1  for  CO-chemistry  by  modifying  the  ReaxFF  angular  terms  and  re¬ 
parameterizing  the  ReaxFF  C/H/O  parameters  against  DFT-data  from  literature. 
High-temperature  ReaxFF  simulation  using  the  new  C/H/O  description  confirm  a 
much  better  agreement  between  intermediates  observed  for  the  CO/H2/O2  system 
between  ReaxFF  and  experiments.  Also,  the  problem  related  to  H-abstraction  by  O2, 
which  is  an  important  initiation  reaction  for  oxidation  of  n-alkanes  have  been  fixed. 
The  newly  developed  force  field  has  retained  its  quality  for  larger  hydrocarbon 
combustions.  However,  the  kinetics  of  alkyl  benzene  pyrolysis  has  not  been 
successfully  captured  by  this  force  field,  which  will  be  subjected  to  further  training. 

Furthermore,  we  have  designed  a  computational  methodology  for  calculating  flame 
speeds  -  providing  a  new  concept  for  comparing  atomistic-scale  simulations  and 
experiment. 
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